COMPUTING CENTRAL VALUES OF L-FUNCTIONS 



1. 

How fast can we compute the value of an L-function at the center of the critical 
strip? 

We will divide this question into two separate questions while also making it 
more precise. Fix an elliptic curve E defined over Q and let L(E, s) be its i-scrics. 
For each fundamental discriminant D let L(E, D, s) be the L-series of the twist Ejj 
of E by the corresponding quadratic character; note that L(E, 1, s) = L(E, s). 

A. How fast can we compute the central value L(E, 1)? 

B. How fast can we compute L(E, D, 1) for D in some interval say a < D < b? 

These questions are obviously related but, as we will argue below, are not iden- 
tical. 

We should perhaps clarify what to compute means. First of all, we know, thanks 
to the work of Wiles and others, that L(E, s) — L(f, s) for some modular form / 
of weight 2; hence, L(E,s), first defined on the half-plane 5ft(s) > 3/2, extends to 
an analytic function on the whole s-plane which satisfies a functional equation as 
s goes to 2 — s. In particular, it makes sense to talk about the value L(E, 1) of 
our L-function at the center of symmetry s = 1. The same reasoning applies to 
L(E,D,s). 

As a first approximation to our question we may simply want to know the real 
number L(E, D, 1) to some precision given in advance; but we can expect something 
better. The Birch-Swinnerton-Dyer conjectures predict a formula of type 

(1) L(E,D,l) = K D m 2 D , 

for some integer nip and kd an explicit easily computable positive constant. (Up 
to the usual fudge factors the conjectures predict that m D , if non-zero, should be 
the order of the Tate-Shafarevich group of Ed ) To compute L(E, D, 1) would then 
mean to calculate mjj exactly. 

In fact, formulas a la Waldspurger have the form {1} with mo the |D|-th coeffi- 
cient of a modular form g of weight 3/2 which is in Shimura correspondance with 
/. The main point of this note is to discuss informally how explicit versions of such 
formulas can be used for problem B above. 

Let us also note the interesting fact that mo, being related to the coefficient of 
a modular form, typically does not have a constant sign. The significance of the 
extra information provided by sgn(mu) remains a tantalizing mystery. 

2. 

There is a standard analytic method to compute L(E, 1), which we now recall. 
If E has conductor N then the associated modular form / has level N and 

f\tO N = ~ £/: 
1 
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where wn is the Fricke involution and e is the sign of the functional equation for 
L(E,s). Concretely, we have 



It follows that 




Now break the integral as J Q + , make the substitution 1 1— > 1/t in the first and 
use the functional equation to obtain 

(This is the classical argument to prove the functional equation of L(E, s) and goes 
back to Riemann who used it for his zeta function.) 
Now plug in s = 1 to get 

(2) ^ i(£i i) = (1+E )^^-W^ 

77 n>l U 

where / has Fourier expansion 

n>l 

Assume for simplicity that gcd(iV, D) = 1. Then the conductor of Ed is ND 2 and 
d2J applied to Ed yields, more generally, 

(3 ) L (E, D. 1) = (1 + e D ) V ( -) ^ e ~ 2 ^ D ^ , 

n>l v 7 

with Ed the sign in the functional equation of L(E, D, s). 

We know that \a n \ grows no more than polynomially with n (a straightforward 
argument gives \a n \ = 0(n)). It follows that for a fixed E and varying D we will 
need to take, very roughly, of the order of 0(|D|) terms in the sum to obtain a decent 
approximation to L(E, D, 1). Assuming the Birch-Swinnerton-Dyer conjectures we 
may use ([3]) to compute m D in ([J) exactly. However, if we know that mu is the D-th 
coefficient of some specific modular form (i.e. we have a formula a la Waldspurger) 
we would get \mo\ but would not be able to recover sgn(mu)- 

Using this method to compute, say, L(E,D,1) for |D| < X would take time 
of the order of 0{X 2 ). We will see below that using formulas of type ([T]) we can 
reduce this to 0(A 3 / 2 ) for at least some fraction of such Z?'s. 

3. 

Before tackling L(E,D,1) let us consider the case of the special value of an 
Eisenstein series of weight 2 (as opposed to a cusp form as we have for L(E, D, 1)). 
What follows is meant only as an illustration of the general case. 

Let the L-function be £((— ) , s— 1)£((— ) , s) with D < the discriminant of an 
imaginary quadratic field K. Its value at s = 1 is essentially h(D) 2 , where h(D) is 
the class number of K, and we find an analogue of ((T|) with h(D) playing the role 
of rriD- There are many excellent algorithms for computing the class number h(D) 
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(see for example [T] chap. 5). Unfortunately, these do not obviously generalize to 
the calculation of mjj. The main reason for this is that the class group of K is 
easy to describe (both its elements and the group operation) in terms of binary 
quadratic forms, whereas its elliptic analogue, the Tate-Shafarevich group of Ejj, 
is notoriously intractable. 

The standard analytic method of the previous section yields the following formula 
(which was known to Lerch, see 2 vol. Ill, p. 171) 



where wd is the number of units in K and cr(n) :— ^2 d < n d is the divisor sum 
function. Again, we need to take, roughly, 0(|Z?|) number of terms in the sum to 
obtain a reasonable approximation of the left hand side. In this case, we in fact 
have an exact formula requiring D terms, namely, Dirichlct's class number formula 



Neither one of these formulas is, however, particularly useful for computing h(D) in 
practice. On the other hand, it may be worth pointing out that similar arguments 
yield the formula [5j vol. Ill, p. 153. 



with the number of necessary steps now reduced to the order of 0(\J\D\). (Anal- 

ogous formulas can be given for D in other congruence classes modulo 8.) 

To make the connection to the general case of computing L(E, D, 1) that we are 
considering we mention two other possible approaches to computing h(D) that do 
generalize. 

(I) The first is to follow Gauss and realize ideal classes of K as classes of primitive, 
positive definite binary quadratic forms of discriminant D. Each class has a unique 
representative Q = (a, b, c) in the standard fundamental domain (what is known as 
a reduced form) and we can simply enumerate these. A straightforward algorithm 
is as follows: run over values of b with b = D mod 2 and < b < y/\D\/3; for each 
b decompose (b 2 — D)/A as ac with < a < c. Add one or two to the total count 
as the case may be if gcd(a, b,c) = 1. 

Though this algorithm also takes time 0(|D|) the constant of proportionality is 
very small making the algorithm quite practical. An important point to notice for 
our purpose, however, is that if we wanted to compute h{D) for < \D\ < X we 
may simply run over all triples a, 6, c of size at most \/a/3 checking the necessary 
conditions on (a, b, c) for it to be a reduced form. In this way we obtain an algorithm 
which will run in time 0(AT 3 / 2 ). 

(II) The second approach is again to follow Gauss but in a different direction. He 
proved that h{D) is related to the number of representations of \D\ as a sum of 
three squares. One precise form of this relation is the following identity (see [3] 






D = 5 mod 8, 
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p.177) 



(G) 



1 

2 



£ q x 2 W+z> = i +12 Y^H 2 (D)qW 




where D runs through all negative discriminants (i.e. D < and D = 0, 1 mod 4), 
and i?2 is a variant of the Hurwitz class number (see [3], page 120). (For us it suffices 
to know that it is related to h(D); for example for D = 5 mod 8 a fundamental 
discriminant we have ^(-D) = h(D).) 

There are sophisticated techniques for computing the coefficients of the left hand 
side, such as convolution which uses the fast Fourier transform to compute products 
of g-series. But even a simple enumeration of the lattice points x 2 +y 2 + z 2 < X, x = 
y = z mod 2 would again take time 0(X i l 2 ). 

The two approaches (I) and (II) are of course related; they amount to counting 
(in an appropriate sense) the number of representations of D by a certain ternary 
quadratic form. In case (I) we count the number of solutions to b 2 — 4ac = D 
up to <SX2(Z)-equivalence; in (II), the number of solutions to \D\ — x 2 + y 2 + z 2 
with x = y = z mod 2. Note the crucial difference that the ternary quadratic form 
involved is indefinite in case (I) and positive definite in case (II). 

A more geometrical point of view is to think that we are dealing with Heegner 
points. In case (I) we may associate to a primitive positive definite binary quadratic 
form Q = (a,b,c) the point zq = (—6 + \/~D)/2a in the upper half plane TL. The 
respective actions of SL2CZ) on forms and 7i are compatible; hence, the class of Q 
determines a unique (Heegner) point in SX2(Z)\7i of discriminant D. 

It is a bit less intuitive how to think of Heegner points in case (II) but this 
was worked out by Gross [3] . The main ingredient is a positive definite quaternion 
algebra B over Q ramified, say, at 00 and a prime N . Pick a maximal order R of 
B and let 1%, . . . ,I n be representatives for the (left) ideal classes of R. Let Ri be 
the right order of Ii for i = 1, • • ■ , n. 

Fix an imaginary quadratic field K of discriminant D. Then we can think of a 
Heegner point of discriminant D (what Gross calls a special point) as an (optimal) 
embedding of the ring of integers Ok into some R i . Eichler has proved that the total 
number of such points, each counted up to conjugation by R*, is (1 — (^)h(D). 
(In fact, the situation is quite analogous to that of case (I) if we take the indefinite 
algebra B = M 2 (Q) and R = M 2 (Z).) 

For example, if N — 2 then the algebra B is the usual Hamilton quaternions and 
we may pick R to be the order discovered by Hurwitz (in standard notation) 



In this case there is only one class of left i?-ideals represented by R itself. Hence a 
Heegner point is an embedding <j> : Ok — * R- 

How do we find such embeddings? The main thing we need is a w G R with 
w 2 = D. Such a quaternion, because D is a scalar, necessarily has trace t(w) = 
and norm n(w) = —D and conversely. Elements of trace in R form a rank 3 lattice 
and hence n(w) = —D is a representation of D by a certain ternary quadratic form 
associated to R. A few congruence conditions are needed to actually produce an 
optimal embedding <p out of w but the upshot is that the problem becomes one 
about representations of — D by ternary quadratic forms. For example, in the case 
N = 2 Eichler's count of embeddings can be completely encoded in the identity 



R = Z + 1A + Zj + Z|(l + i + j + k). 
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([6]); the presence of the factor 12 in that formula is due to the fact that this is the 
order of R* / ± 1. More details on this setup are given below in §4 (II). 

4. 

We now return to the main case of computing L(E, D, 1) and describe analogues 
of cases (I) and (II) of the previous section. These analogues are the remarkable 
results of Gross and Zagier. 

(I) Let us assume for simplicity that E has conductor a prime N, sign of the 
functional equation equal to —1, L'(E, 1) 7^ 0, and E(Q) = (Pq). If / is the weight 
2 eigenform associated to E then we get a map 

$: X (N) — C/L 

z 1 ► 2ni J ioo f(u) du 

where X (N) is the modular curve of level N and L c C is a certain lattice of 
periods of /. It is known that C/L = E'(C) for some elliptic curve E' /Q isogenous 
to E. Since the L-function is unchanged by isogenies we may assume without loss 
of generality that E' = E. 

Let K be an imaginary quadratic field of discriminant D < —4 in which N 
splits. Choose 6* £ Z such that 6* = D mod N; this is possible by the assumption 
that N splits in K. Note also that N does not divide 6*. We want to consider 
Heegner points on Xo(N) of discriminant D. To define them concretely choose 
representatives Q = (a, b, c) of the h(D) classes of binary quadratic forms with 
N I a and 6=6* mod N. (For example, start with representatives (a, b, c) with 
gcd(a, N) — 1 and compose them with the fixed form (N, &*, (6» - D)/2N).) 

Then z Q := (-6 + VD)/2/a) e X (iV) is well defined and P D := Eg $ 0q) € 
E(K). Moreover, complex conjugation fixes Pq, by the assumption on the sign of 
the functional equation. Hence Pd actually is in E(Q) (and is independent of the 
choice of 6*). 

One consequence of the results of Gross-Zagier is the following [T5] , [3] , [5] . By 
our assumption on E(Q) we have Pd = moPo for some raufiZ and hence 



(8) L(E,D,1) = KD m 2 D ; 

where kd is an explicit easily computable positive constant; i.e. we have a formula 
of type (HJ). 

Usually one regards the Gross-Zagier formula as a way to compute a rational 
point Pd on E whose height is given in terms of L(E, D, 1)L'(E, 1) and hence 
obtaining, when this value does not vanish, a confirmation of the predictions of the 
Birch-Swinnerton-Dyer conjecture. Here, instead, we are taking the point of view 
that the points of E(Q) are known and use the Gross-Zagier formula as a means 
to computing L(E, D, 1). 

To calculate mjj in practice it is better to work on the E(C) = C/L model of 
E rather than, say, a Weierstrass equation. Let zq € C represent, modulo L, the 
point Pq £ E(Q). We first compute an approximation to 



EE? 



z D := > > _ e 2 «» z «. 



Then we solve the linear equation below for integers n\ and 712 

zq = mDZo + mtjJi + n 2 uJ2, 
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where uj\,uj2 are a basis for L. (In fact, multiplying by 2 if necessary we may 
assume that u)\ € R and u>i £ iM. and hence by taking real parts solve only a three 
term equation instead.) 

The result is a practical and reasonably efficient algorithm for computing tud- 
The number wb is the D-th Fourier coefficient of a weight 3/2 modular form g of 
level AN which is in Shimura correspondence with /. It is interesting that we can 
compute the Fourier coefficients of g directly without any knowledge of the whole 
vector space of modular forms in which g lies; though we do, of course, start by 
knowing / itself. (We have only described the calculation for certain D's but there 
is analogous way to get all coefficients.) 

Together with my student Ariel Pacetti we implemented the above algorithm in 
GP. The corresponding routines can be found at 

\protect\vrule widtriOpt\protect\href {http : //www. ma . utexas . edu/users/villegas/cnt/}{http : //wwi. 

under Heegner points. 

Here is a sample example. Let E be the curve y 2 + y — x 3 — x of conductor 
N — 37 (this is the elliptic curve over Q of positive rank with smallest conductor) . 
This case was described in detail in [TS]. It is known that E(Q) — ((0, 0)). 

? e=ellinit( [0,0, 1,-1,0] ) ; anvec=ellan(e,5000) ; 

? for(d=5,100, if (isf undamental (-d) && kronecker (-d,37)==l , 

print (-d, " " , ellheegnermult (e , -d, [0,0] ,0,anvec) [1] ))) 

-7 -11 -40 -47 -67 -71 -83 -84 -95 
1-1-21-6-11 10 

The first row is D, the second raj (for typographical reasons we transposed the 
actual GP output). These values agree, fortunately, with Zagier's [15] formula (28) 
up to a global negative sign. 

In our implementation at least the algorithm is not that well suited for computing 
L(E, D, 1) for all D < and \D\ < X for very large X; for this, it would be better 
to adapt (see §5) the ideas of (II) below but these have not been fully implemented 
as yet. 

(II) Let B over Q be the (unique up to isomorphism) positive definite quaternion 
algebra ramified at oo and a prime N . Pick a maximal order R of B and let 
I\, . . . , I n be representatives for the (left) ideal classes of R. Let Ri be the right 
order of li for i = 1, • • ■ , n. The class number n of i?, in contrast with h(D), has a 
simple formula and is roughly of size N/12. 

For example, if N = 3 mod 4 we can describe B as the algebra over Q with 
generators i, j such that i 2 — ~l,j 2 = —N and ij = —ji. Also in this case we can 
take R = Z + Zi + Z±(l + j) + H\(\ + j). 

There are various ways to compute representatives I\, ■ ■ ■ , /„ of the ideal classes 
(for algorithms for quaternion algebras see [H]). If N = 3 mod 4 there is an 
algorithm which is completely analogous to that of Gauss §3 (I) for binary quadratic 
forms. It exploits the fact that our choice of R has an embedding of Z[i] and hence 
allows us to view i?-left ideals as rank 2 modules over Z[i]; then classes of R- 
ideals correspond to classes of positive definite binary Hermitian forms over Z[i] of 
discriminant —N. Instead of H. we now need to work on hyperbolic 3-space where, 
as it turns out, the action of SX2(Z[i]) has a very simple fundamental domain. This 
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yields an algorithm which is almost verbatim that of Gauss for binary forms over 
Z. Details can be found in [12] , 

For example, if N — 11 then there are two classes of positive definite binary 
Hermitian forms of discriminant —11 over Z[i]; namely, (1,1,3) and (2,1 + 2z, 2) 
corresponding to the two ideals 

Jo := R = Z + Zi + Z\{l+i)+Zi\{l+j) 

and 

Ji := 2Z + Z2i + Z±(l + 2% + j) + Zi\{\ + 2i + j) 

representing the n — 2 classes of left i?-ideals. 

Let Vq be the Q vector space of functions on the set {Io, ...,/„}. For each 
m G Z>o there is an operator B(m) acting on Vq, the Brandt matrix of order to, 
which encodes the number of representations of to by certain quaternary quadratic 
forms (see [3] (1.4)). Let B be the algebra generated over Z by all the B(m); it is 
commutative and B®zQ is semisimple. 

On the other hand, we have the space Mc of modular form of weight 2 on T (N) 
(known to be of dimension n) and the Heckc operators T m acting on Mc- Let T 
be the algebra spanned by the T m over Z ; like B it is commutative and T ®z Q is 
semisimple. This algebra preserves the Q vector space Mq C Mc of dimension n 
consisting of those modular forms in Mc with Fourier coefficients in Q. 

These two setups are closely related and indeed we have a special case of the 
Jacquet-Langlands correspondence. Eichler proved that T m and B(m) have the 
same trace for all to £ N. Hence, by semisimplicity of the algebras the map T m i— > 
B(m) induces a ring isomorphism T ~ B. It follows that eigenspaces of Vq and 
Mq, under the action of B and T respectively, correspond to each other. Since we 
also have multiplicity one these eigenspaces are one-dimensional. 

In conclusion, given / = 53«>o a nQ n £ Mc an eigenform for all Hecke operators 
T m (so that T m f = a m f) there is an e/ £ Vq ®i K unique up to scalars such that 
B(m)ef = a m ef. (Here K denotes the field Q(ao, oi, . . .) generated by the Fourier 
coefficients of /.) 

In fact, this correspondence gives an efficient way to compute Fourier coefficients 
of eigenforms in Mc (see |llj). An implementation of the corresponding algorithms 
can be found in the above mentioned website (under qalgmodf orms). Here is a 
sample GP session. 

? R=qsetprime (11) ; 
? brandt(R,2)~ 
[1 3] 
[2 0] 

? brandt(R,3)~ 
[2 3] 



[2 1] 
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The first line defines R as a maximal order in the algebra ramified at 1 1 and oo ; the 
others compute the corresponding Brandt matrices. We find that these matrices 
have two eigenvectors: = (1/2,1/3) and e/ = (—1,1) corresponding to an 
Eisenstein series and a cusp form, respectively. 

The above implementation is intended for small to medium scale computations. 
For large scale computations one should use the graph method ideas of Mestre and 
Oesterle jS], which exploit the sparse nature of the Brandt matrices. 

Now following Gross we show how to associate a modular form of weight 3/2 to 
an eigenvector e/. Let Ri be the right order of Ii and let Li c Ri be the ternary 
lattice defined by 

Li : w € Ri, t(w) = 0) mod 2Ri. 

Let gi be the corresponding theta series 

0i(r) := \ J2 9 n(w) , q = e 2 ™. 
weLi 

Gross [3] prop. 12.9 describes precisely how the D-th coefficient a,(D) of gi relates 
to the optimal embeddings of imaginary quadratic orders of Q{yD) into Ri. 

These theta series are modular forms of weight 3/2 and level AN and, in fact, 
belong to a certain subspace U defined by Kohnen. This subspace is determined by 
the condition that the coefficient of q of a form should be zero unless D := —d is a 
discriminant, i.e., D = 0, 1 mod 4, and (-^) 7^ 1. The weight 3/2 Hecke operators 
T m 2 preserve U. 

Define 

9 : = £ e /(*) 9i = J2 mD 9 l ° l e U - 

i D 

This form is identically zero if the sign in the functional equation of / is — 1. If g is 
non-zero it is a modular form in Shimura correspondence with /; i.e., T m ig — a rn g, 
where T m f = a m f. Moreover, we have the Waldspurger formula [3] 13.5 

(9) L{f,l)L{f® XD ,l) = K f ^=m 2 D , 

where D is a fundamental discriminant with (^) 7^ 1, xd is the associated qua- 
dratic character, Kf > is a constant depending only on / and Sd := 2 if N \ D 
and Sd '■= 1 otherwise. 

Finally, let E/Q be an elliptic curve of prime conductor N and sign +1 in its 
functional equation. Let / and g be the corresponding modular forms of weight 2 
and 3/2 respectively as above. Then if L(f, 1) 7^ we obtain from © a formula 
of type ([1]) with itid the Fourier coefficient of g. As in §3 (II) to compute tud for 
\D\ < X we could run through all w S Li with n(w) < X whose total number 
is 0(X 3 ' 2 ). Again various computational techniques could also be used to speed 
up the calculation of tod. Note that in any case all computations are done with 
integer arithmetic. 

Tables of m,D 's for several curves and the routines to compute them can be found 
at G. Tornana's website 

\protect\vrule widthOpt\protect\href {http : //www . ma . utexas . edu/users/tornar ia/ cnt/}{http : //wi 
among other goodies (an interactive version of Cremona's tables of elliptic curves 
and an interactive table of ternary quadratic forms). 
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■5. 

We conclude with some remarks about the general situation. 

1. It follows from © that if L(f, 1) = then the form g vanishes identically. In 
this case we naturally need to do something else. 

In [7] we work out an extension of Gross's work introducing an auxiliary prime 
I; the theta series gi, for example, are modified by introducing an appropriate 
weight function. The complexity of algorithms only increase by a factor essentially 
proportional to 

2. If the level N is not prime but square- free the situation is not too different 
from the one described above. The downside is that L(E, D, 1) can be computed 
this way only for a certain fraction of D's (determined by local conditions). One 
needs to consider a quaternion algebra B ramified at oo and at primes I | N for 
which the Atkin-Lehncr involution acts as f\ Wl = —f and an Eichler order in B of 
level the product of the remaining primes factors of N. 

3. If the level is not square-free things become quite a bit more complicated; 
for example, the algebra B of Brandt matrices typically does not act with mul- 
tiplicity one and some modular forms are simply missing. The arithmetic of the 
corresponding orders, which are no longer Eichler orders in general, also becomes 
more involved and, moreover, one needs to consider two types of orders: one for 
the weight 2 side and another for the weight 3/2 side; see [9], [10], [14] some work 
on this case. 

4. To compute twists L{f®xu 1) by rea >l quadratic fields Q(Vl) one may consider 
a twist fjj := f <S> xd by an auxiliary imaginary quadratic field Q(vD) and find 
a formula of type ([I]) for L(fo <8> xdiA)- The form fr> typically does not have 
square- free level so several corresponding difficulties ensue, see [14]. 

5. Forms of higher weight can also be handled using quaternion algebras by 
introducing harmonic polynomials as weight functions for the theta functions (both 
for the ideals Ii corresponding to forms of weight 2 + 2r and for the ternary lattices 
Li corresponding to forms of weight 3/2 + r) see [6], [13]. 
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